Thermal convection in fluidized granular systems* 
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Thermal convection is observed in molecular dynamic simulation of a fluidized granular system of 
nearly elastic hard disks moving under gravity, inside a rectangular box. Boundaries introduce no 
shearing or time dependence, but the energy injection comes from a slip (shear-free) thermalizing 
base. The top wall is perfectly elastic and lateral boundaries are either elastic or periodic. The 
observed convection comes from the effect of gravity and the spontaneous granular temperature 
gradient that the system dynamically develops. 
PACS numbers: 45.70.Mg, 47.20.Bp, 47.20.Te, 81.05. Rm 



In the study of granular systems, convection has at- 
tracted particular attention. Many of the experimental 
and simulation studies of granular convection focus their 
attention on the effects of a vibrating base ■ Par- 

ticular care has been taken to determine the different 
patterns that appear depending on the amplitude and 
frequency of the vibrating base, the form of the interac- 
tion force between the grains, the effects caused by the 
roughness and inclination of the walls and the role of 
voids in this vibrating-base-plus-wall convection. Some 
theories based on hydrodynamic continuum equations for 
this type of convection have been developed |K|-|T2]] . Also 
the effect of the internal shear bands in th e sy stem has 
been pointed out as a source of convection [13|] . 

The purpose of this letter is to point out the existence 
of a convective regime in granular systems when there is 
no vibrating base, or any effect caused by the roughness 
or the angle of the walls, or any shearing effect caused 
by the walls themselves. The source of this convection 
stems from gravity and the dissipative nature of the gran- 
ular collisions, and some hint about its existence was al- 
ready mentioned in In the present model system, each 
wall is represented by a shear-free and time-independent 
boundary condition. 

More precisely the bottom wall gives a stochastic nor- 
mal component to the velocity of each particle that col- 
lides with it. The tangential velocity is unchanged, thus 
a shear-free thermal boundary condition is imposed. The 
top wall is perfectly elastic while the lateral walls are ei- 
ther perfectly elastic or they correspond to a periodic 
boundary. Particles are subjected to the acceleration of 
gravity g. 

When energy is pumped into a granular system 
through a boundary — the bottom wall in the present 
case — a temperature gradient develops spontaneously. 
Since the particle-particle collisions are dissipative the 



system is hotter close to the injecting boundary and 
colder away from it. There is then an energy flux from 
the base upwards which is dissipated in the bulk through 
collisions. If the system is almost perfectly elastic it re- 
mains macroscopically static. 

Performing molecular dynamic simulations of inelastic 
hard disk systems in the presence of gravity, a transition 
from a hydrostatic (purely conductive regime) to a con- 
vective regime is observed by increasing the dissipative 
parameter present in the particle-particle collision rule. 

In standard fluids the onset of thermal convection is 
roughly determined by the ratio between the character- 
istic times of the processes against convection (viscous 
and thermal diffusion) and favorable to convection (buoy- 
ancy), through which the externally imposed tempera- 
ture gradient plays a role H]. The idea that underlies 
this letter is that in a granular system there are similar 
mechanisms which lead to a transition from a conductive 
to a convective regime. 

There is an important difference, though. The tem- 
perature gradient is not externally imposed but rather it 
is dynamically created by the system itself. This means 
that whenever energy is injected into the system a tem- 
perature gradient develops (due to the dissipative nature 
of the particle-particle collisions) and not three but four 
ingredients compete inhibiting or fostering the appear- 
ance of convection: viscous and thermal diffusion, buoy- 
ancy and dissipation. 

This letter reports the results about convection ap- 
pearing in a 2D system of N hard disks with mass m = 1 
and diameter a — 1, which collide inclastically with the 
rule 



"12 



= -(1 - 2q) (vi2 ■ n) , vu'-t=vu-i, (1) 



where V12 = c± — c-i is the relative velocity between the 
colliding particles, the primed and unprimed variables re- 
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fer to the post and precollisional velocities, ft and i are 
the unit vectors normal and tangential to the contact 
plane, and q is the dissipative coefficient q = (1 — r)/2, 
r being the normal restitution coefficient. Only trans- 
lational degrees of freedom are present. The simulation 
results reported below come from our event-driven molec- 
ular dynamic simulations |l5[ , and the careful measuring 
routines developed in p6| . 

As already mentioned, the system is maintained in a 
fluidized state by the injection of energy from a thermal 
(stochastic) base (y = 0) at temperature T(0) = Tbasc = 
1 in energy units, while the top boundary is a perfectly 
elastic wall at y — L. Gravity enters the problem through 
the Froude number, 



Fr 
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Although there is a wide range of parameters for which 
convection occurs, all the simulations being reported cor- 
respond to N — 2300 particles, fraction of occupied area 
Pa = 0.18, and relative gravity strength Fr = 0.1. The 
only parameter that varies from one simulation to an- 
other is q, even though for convenience the product qN 
is used in what follows. 

a. Convection inside a box. Consider first the granu- 
lar system inside a squared box with a thermal base, and 
perfectly elastic upper and lateral walls. As the elastic 
walls do not produce any direct shear on this system and 
the bottom base is stochastic but preserving the hori- 
zontal component of the velocity, the boundaries intro- 
duce neither spatial nor temporal macroscopic correla- 
tions. It can be said that none of the usual conditions 
under which convection has been studied in granular sys- 
tems are present, nevertheless convection does appear, as 
observed in Fig. m. 
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FIG. 1. Averaged velocity field for different values 
of qN. At the left qN = 6 and at the right qN = 34. 

The hydrodynamic stationary solution for the associ- 
ated conservative system [q = 0) is simply a constant 



temperature system with no heat flux and density de- 
creasing with height. Including a small amount of energy 
dissipation in each collision, one and even multiple rolls 
may develop into the system. Hence the convection we 
observe is due solely to gravity and the dissipative nature 
of the collisions between particles. 

To understand the origin of this convection, we have 
plotted (see Fig.^J) the dependence on the dissipative pa- 
rameter qN of the observed granular temperature dif- 
ference between the bottom (y — 0) and the top of the 
system (y = L), A = T(0) - T(L). For small values 
of qN this difference increases with increasing qN. This 
situation corresponds to the conductive regime, in which 
A can be written as an expansion on the small param- 
eter qNpA [0- This difference A reaches a maximum, 
at about qN w 4, and from then on it decreases, prov- 
ing that a mechanism favoring energy transport from the 
base upwards appears at this value. This fact is also cor- 
roborated by the amount of energy per unit time (call it 
heat flux) Q(0), entering the system through the base, 
plotted in the same figure in arbitrary units. It is seen 
that the heat flux is steeper about the same qN for which 
A reaches its maximum: more energy per unit time is re- 
quired to keep stable the temperature at the base. 

From these observations we can conclude that there 
exists a threshold qN value from which the convection is 
triggered. Furthermore, this also supports the idea that 
convection starts when a critical value of the temperature 
difference between the bottom and the top of the system 
is reached. 
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FIG. 2. Differences between the bottom and top 
temperatures (solid circles) and rescaled heat flux (open 
circles) versus qN 

A way to detect and quantify this transition is by mea- 
suring mass circulation in the system. This is imple- 
mented by calculating the sum of integrals of the veloc- 
ity field along many concentric paths centered about the 
geometric center of the box: $ = ^ J v-dl. This observ- 
able <i> will be negligible if there is no convection and it 
will be distinctly nonzero (positive or negative) if there is 
one (anti) clockwise convective roll. Observed values of i> 
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are plotted in Fig. ^ clearly showing a supercritical tran- 
sition at qN ~ 4 from the conductive to the convective 
regime with one convection roll. Due to the symmetry of 
the problem, rolls with both signatures are equally prob- 
able and they appear in our simulations, depending only 
on the initial condition. 

It can also be seen that from qN w 34 up there is 
a coexistence of regimes with zero and non zero circula- 
tion which corresponds to the competition of one and two 
convective rolls, this is, a subcritical transition from the 
one-roll to the two-rolls regime. As the dissipation coef- 
ficient q continues increasing, other transitions to multi- 
rolls patterns can be observed with difficulty but since 
the system starts getting denser all convective movement 
eventually disappears. 
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FIG. 3. Mass circulation $ measured in simula- 
tions (points). The dashed lines correspond to the curve 
$ = ±0.04VgA r - 3.8. 

Although no theory has been presented yet, it seems 
that the appearance of the two and multi-rolls regimes 
could be due to a change in the effective aspect ratio of 
the system: as dissipation increases, there are regions 
where density rises considerably, lowering the average 
height occupied by the system. 

It is worth mentioning that when two rolls were ob- 
served they always appeared as shown in Fig. |l|, namely, 
the fluid goes up in the middle of the box and comes 
down along the walls. This privileged signature seems to 
have its origin in the local increase of density that walls 
cause. Higher density implies more collisions and there- 
fore more dissipation, hence lower granular temperature: 
the system is heavier near the walls. 

b. Convection with horizontal periodic boundary con- 
ditions. Any effect that the elastic lateral walls could 
have on the onset of convection in the previous case is 
discarded when periodic lateral boundary conditions are 
imposed on the system. 

The container is a periodic channel, and in this case 
a transition to a convective regime is found again, al- 



though, due to the absence of lateral boundaries, the 
convective rolls appearing in the system travel now along 
the channel. This was observed even though the simula- 
tions were carefully initialized with zero total horizontal 
momentum P x , namely with zero horizontal mass flux. 
Since the boundary conditions do not change the horizon- 
tal component of the velocities, P x remains zero during 
the evolution, as was confirmed in the simulation. 

In order to detect this pattern, we performed time av- 
erages of the mass flux field. Because of the roll move- 
ment, this averaging time must be larger than the mi- 
croscopic time and smaller than the time needed for the 
roll to travel a significant distance. We chose this time 
to be much smaller than the ideal gas thermal diffusion 
time which, in our units, is of order N^y^^T\^ c , but 
large enough to contain multiple particle-particle colli- 
sions. The observed rolls persisted for times longer than 
the macroscopic time, resulting in an hydrodynamic pat- 
tern. 
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FIG. 4. Mass flux field averaged in circles of 250 
collisions per particle. The figures a), b), c) and d) cor- 
respond to cycles 100,112,124 and 136 respectively. A 
big roll can be observed moving to the left side of the 
system while a small roll appears varying its size. 

An example of what is happening in the system can 
be observed in figure ^. This figure is a plot of the aver- 
aged mass flux field at four different stages of the simula- 
tion. Due to the periodic lateral boundaries, the solution 
should be a two-rolls pattern, (or any even number of 
rolls) but the aspect ratio forced on the system would 
imply rolls with a width about half their height, which 
makes them unstable. The system was most of the time 
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observed to have one large roughly circular roll accom- 
panied by a smaller one. 

The reason for this movement seems to be found in 
vortex dynamics. It is well known that a vortex near a 
fixed isolating wall behaves as if it were in front of a twin 
vortex which ensures the condition of null hydrodynamic 
velocity v y — at the top elastic wall This twin 

vortex would induce a movement parallel to the wall and 
a sense determined by the sign of the circulation of the 
original vortex. 



In conclusion it can be stated that bidimensional gran- 
ular systems present convective regimes when there is 
gravity and a thermal base even though no shearing is in- 
troduced through the boundary conditions. It has been 
argued that such convection owes its existence only to 
gravity and the dissipative nature of the particle-particle 
collisions. If the system approaches the clastic limit such 
convection disappears as there is no spontaneous creation 
of a vertical temperature gradient. This gradient plays 
the key role in this kind of convection. 




FIG. 5. Total velocity correlation Ctot vs. qN. The 
transition from a conductive to a convective regime 
is observed at qN ~ 18. The data have some dis- 
persion near the transition due to the finite size of 
the system. The dashed line corresponds to the curve 
Ctot = 0.00137^ - 18. 

In the present case, to detect the onset of convection, 
it was not practical to measure the circulation $ as when 
the system was between hard lateral walls. Instead the 
transition was detected measuring a space velocity corre- 
lation. The system is tiled with cells and a hydro- 
dynamic velocity correlation is defined by 



C{i,j) = \Y,^h3)-v{i',3'). 



(2) 
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where the cells («', j') refer to the eight first neighbors 
of the cell. This observable has the advantage 

of being insensitive to the displacements of the convec- 
tive pattern. The total correlation is defined as Ctot = 
l/A^us^T C(i,j). It measures how similar, on the av- 
erage, are the velocities in neighboring cells. When there 
is no convective pattern in the system, then the time 
averaged value of Ctot is nearly zero, while as soon as a 
convective current develops, Ctot takes distinctly positive 
value. 

Figure |)| shows the evolution of C to t with qN. The 
transition from conductive to convective regime is clearly 
seen. It takes place roughly at qN =18. 
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